{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8bad2592-8d90-4a80-9123-42f51a89b346",
   "metadata": {},
   "source": [
    "# [Daily Dose of Data Science](http://join.dailydoseofds.com/)\n",
    "\n",
    "Code accompaning the newsletter issue: [Implementing a Siamese Network](https://blog.dailydoseofds.com/p/implementing-a-siamese-network)\n",
    "\n",
    "Author: Avi Chawla"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "248d8978-b082-432f-ad47-ee784dc6f255",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "45272686-9178-4dcd-a9c3-f0bb826ae82a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import numpy as np\n",
    "import random\n",
    "import torchvision.transforms as transforms\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from torch.utils.data import DataLoader, Dataset\n",
    "from torchvision.datasets import MNIST\n",
    "from torch import optim"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a483933c-0fe6-4b02-87de-b3ca28dcc2a7",
   "metadata": {},
   "source": [
    "## Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "e924e567-a393-4df5-bb5c-97bd55701af6",
   "metadata": {},
   "outputs": [],
   "source": [
    "class SiameseDataset(Dataset):\n",
    "    def __init__(self, data, transform=None):\n",
    "        self.data = data    \n",
    "        self.transform = transform\n",
    "            \n",
    "    def __getitem__(self, index):\n",
    "        imgA, labelA = self.data[index]\n",
    "        \n",
    "        same_class_flag = random.randint(0, 1)\n",
    "        \n",
    "        if same_class_flag:\n",
    "            labelB = -1\n",
    "            \n",
    "            while labelB != labelA:\n",
    "                imgB, labelB = random.choice(self.data)\n",
    "                \n",
    "        else:\n",
    "            labelB = labelA\n",
    "            \n",
    "            while labelB == labelA:\n",
    "                imgB, labelB = random.choice(self.data)\n",
    "\n",
    "        if self.transform:\n",
    "            imgA = self.transform(imgA)\n",
    "            imgB = self.transform(imgB)\n",
    "            \n",
    "        return imgA, imgB, torch.tensor([(labelA != labelB)], dtype=torch.float32)\n",
    "\n",
    "    def __len__(self):\n",
    "        return len(self.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "79241ce5-8ed2-46f5-8333-295c02cbbb89",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mnist_train = MNIST(root='./data', train=True, download=True)\n",
    "mnist_test = MNIST(root='./data', train=False, download=True)\n",
    "\n",
    "transform = transforms.Compose([transforms.ToTensor()])\n",
    "\n",
    "siamese_train = SiameseDataset(mnist_train, transform)\n",
    "siamese_test = SiameseDataset(mnist_test, transform)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a05693f-3516-4cba-a74e-864e0667f7da",
   "metadata": {},
   "source": [
    "## Network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "c2f04b78-33de-4020-a87d-d651d320af07",
   "metadata": {},
   "outputs": [],
   "source": [
    "class SiameseNetwork(nn.Module):\n",
    "    def __init__(self):\n",
    "        super(SiameseNetwork, self).__init__()\n",
    "\n",
    "        self.cnn = nn.Sequential(\n",
    "            nn.Conv2d(1, 64, kernel_size=5, stride=1, padding=2),\n",
    "            nn.ReLU(inplace=True),\n",
    "            nn.MaxPool2d(2, stride=2),\n",
    "\n",
    "            nn.Conv2d(64, 128, kernel_size=5, stride=1, padding=2),\n",
    "            nn.ReLU(inplace=True),\n",
    "            nn.MaxPool2d(2, stride=2),\n",
    "\n",
    "            nn.Conv2d(128, 256, kernel_size=3, stride=1, padding=1),\n",
    "            nn.ReLU(inplace=True),\n",
    "            nn.MaxPool2d(2, stride=2)\n",
    "        )\n",
    "\n",
    "        self.fc = nn.Sequential(\n",
    "            nn.Linear(256 * 3 * 3, 1024),\n",
    "            nn.ReLU(inplace=True),\n",
    "\n",
    "            nn.Linear(1024, 256),\n",
    "            nn.ReLU(inplace=True),\n",
    "\n",
    "            nn.Linear(256, 2)\n",
    "        )\n",
    "\n",
    "    def forward_once(self, x):\n",
    "        output = self.cnn(x)\n",
    "        output = output.view(output.size()[0], -1)\n",
    "        output = self.fc(output)\n",
    "        return output\n",
    "\n",
    "    def forward(self, inputA, inputB):\n",
    "        outputA = self.forward_once(inputA)\n",
    "        outputB = self.forward_once(inputB)\n",
    "        return outputA, outputB"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1663bd05-585a-4b6f-9e0e-64ab74e8062c",
   "metadata": {},
   "source": [
    "## Loss function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "911be2cb-e203-4e0e-afae-5fc9a6747cb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "class ContrastiveLoss(torch.nn.Module):\n",
    "    \n",
    "    def __init__(self, margin=2.0):\n",
    "        super(ContrastiveLoss, self).__init__()\n",
    "        self.margin = margin\n",
    "\n",
    "    def forward(self, outputA, outputB, label):\n",
    "        euclidean_distance = F.pairwise_distance(output1, output2, keepdim = True)\n",
    "\n",
    "        same_class_loss = (1-label) * (euclidean_distance**2)\n",
    "        diff_class_loss = (label) * (torch.clamp(self.margin - euclidean_distance, min=0.0)**2)\n",
    "    \n",
    "        return torch.mean(same_class_loss + diff_class_loss)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67f53c12-3ba8-4521-807c-3cb3d9ca8b58",
   "metadata": {},
   "source": [
    "## Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "833a9119-0e73-4147-9ab5-590147393bb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_dataloader = DataLoader(siamese_train,\n",
    "                        shuffle=True,\n",
    "                        num_workers=8,\n",
    "                        batch_size=64)\n",
    "\n",
    "net = SiameseNetwork().cuda()\n",
    "criterion = ContrastiveLoss()\n",
    "optimizer = optim.Adam(net.parameters(), lr = 0.001 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "4d28c979-c3c0-40c7-98e4-18d8348e6167",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 0; Loss 294.2397748604417\n",
      "Epoch 1; Loss 102.84626789577305\n",
      "Epoch 2; Loss 63.36508319573477\n",
      "Epoch 3; Loss 44.61645963555202\n",
      "Epoch 4; Loss 32.77574067004025\n"
     ]
    }
   ],
   "source": [
    "for epoch in range(5):\n",
    "    total_loss = 0\n",
    "    \n",
    "    for imgA, imgB, label in train_dataloader:\n",
    "\n",
    "        imgA, imgB, label = imgA.cuda(), imgB.cuda(), label.cuda()\n",
    "        optimizer.zero_grad()\n",
    "        outputA, outputB = net(imgA, imgB)\n",
    "        loss_contrastive = criterion(outputA, outputB, label)\n",
    "        loss_contrastive.backward()\n",
    "\n",
    "        total_loss += loss_contrastive.item()\n",
    "        optimizer.step()\n",
    "\n",
    "    print(f\"Epoch {epoch}; Loss {total_loss}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80dbf5b2-05cf-40b2-a5ed-edf8aa8878e9",
   "metadata": {},
   "source": [
    "## Plot outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "ca54774f-7c0f-48fa-9149-9a07c7dacf0c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9199\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWUAAADOCAYAAADmKRMiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAbB0lEQVR4nO3deVxU5f4H8A+oDMgygMSwiaCYXjUoXEC9dRVcckGvpWJqmlr5awFKBS0tl1vuWoq5peaWWuaCpveae1ZKuZBLZUmQXZNB9LImi3B+f3gb7/QcdQYG5hn4vF+veb3kwzlzngNfvx7PM+ccu+LSWwqIiEgK9tYeABER3cGmTEQkETZlIiKJsCkTEUmETZmISCJsykREEmFTJiKSCJsyEZFE2JSJiCTCpkxEJBE2ZTOsX78Ojg71cerUSWsPpVqtXLEcQ4fEIqRZMBwd6uPZMaOtPSQyUV2o0V9//RVv/WMG/topEj7eXvD31aF7tygcPHjA2kOzCDZlEsyfPw9HjhxGq1atUL9+fWsPh8jIp7t3YcH8eWjaLATTps/Aa69PRmFBIfr0ehzr1q219vCqjH/jSLD/wCEEBgbCzs4OjTy01h4OkZG//a0LfkrPgJeXlyF77vmx6NCuLWZMn4aRI5+x3uAsgEfKVfTsmNFo5KHF5cuXMeDv/dDIQ4umQYFYvmwpAOD8uXPo2aMbPN3d0DykKbZs3my0/o0bNzBpYiLaPvIwGnlo8UAjD/SL6YOz334rbOuXX37BkwP+Dk93NzT290XihHHY/9k+ODrUx9GjR4yW/frrVMT07Q1vL094aF3RLborvvrqS5P2qUmTJrCzs6vcD4SkU9tqtFXr1kYNGQA0Gg0e7/U4rvz73ygoKDDzJyQXNmULKC8vR/+YvggIaIyZs2YjsEkTvJIQj/Xr16FfTB+0DW+Lt2fOgquLK8aMfgYZGRmGdTMyfsauXbvQq3dvzJk3H6+OG48L58+je7co/Pbbb4blioqK8HjP7jh06CBefOllTJz0Gk4cP47Jr78ujOfw4UPoFtUVBfn5mDzlDcz4x1vIy83D4z2645tvvq6RnwnJpS7UaFaWHg0bNkTDhg0rtb40iktvKXyZ9lq5arUCQPny+AlDNvzpEQoAZcY/3jJkWdk5ipOTk2JnZ6ds2LjJkJ89d0EBoEye8oYhyysoUn4vLjXazg8/XlI0Go3y5tRphmzO3HkKAGXrJ9sNWW5+odKiRUsFgLJv/wGluPSWcrOkTAkJaa5079FDuVlSZlj2P3kFSlBwsBLdrZtZ++zs7KwMf3qE1X/2fLFG7/W68N0PiqOjozJ02HCr/w6q+uKRsoWMGj3G8Gd3d3c8+GALODs7Y+CgQYb8wRYt4O7ubnQUotFoYG9/+9dQXl6O69evw8XFBQ8+2AJpZ84Ylvvss33w9/dH35gYQ+bo6IjRY+5sFwC+TUvDpUs/ITb2KVy/fh05OTnIyclBUVERunaNwhfHjqGiosLi+0/yq601+vvvv2PoU0Pg5OSEt96eafoPRFKc6LMAR0dHPPDAA0aZm9YN/v4BwrlZN60Wuf/5j+HriooKLElejBUrliMzIwPl5eWG73k28jT8+fIvlxHctKnwfs2ahRh9fenSJQDAs2NG3XW8eXl58PDwMHHvqDaorTVaXl6Op4cNxffff4eU3Z/Cz8/vvuvIjk3ZAurVq2dWruDOE7jmzJ6F6dOmYuQzozB16nR4eHrA3t4eiePHV+qI9o91Zs2eg9CwMNVlXFxczH5fsm21tUZf+L+x2Lt3D9au24CuXaPMHouM2JStbMf27fhbly5YsfJ9ozwvLxeNvBoZvg5sEogfvv8eiqIYHYmkp18yWq9ps6YAAFc3N0RHd6vGkVNdIWuNvjYpCevXrcX8BQsRO2RIpd9HNjynbGX16tWDohg/u3bbJ5/gypUrRln37j1w5coVfLp7tyErLi7GmtWrjZYLD2+Lps2a4d2FC1FYWChs79q1axYcPdUFMtbowgXz8c7ChUiaOAkvx8WbszvS45GylfXq3Rsz334Lzz07BpEdO+LC+fPYsnkTgps2NVru2eeex7JlSzHi6WF46eU4+Pr6YsvmTXB0dAQAw5GJvb09li1fgf4xffHIw6EYMWIk/Pz98duVKzh69AjcXN2wfWfKPce059PdOHv2LACgrKwM58+dw6yZbwMA+vaNwUOhoZb+MZDEZKvRlJ078fprkxAS0hwtW/4Fmz780Oj70d26QafTWfinUHPYlK1s4qTX8PvvRfhoyxZ8svVjPPzII9iRsgtTJht/ttPFxQX/2rcf415NwHtLkuHi4oJhw4cjMrIThsQOMhQ+cPuKp6Off4FZM9/G8mVLUVhYCJ2PDzq074Axzz133zHt2LEDGzesN3ydlnYGaWm3Z9n9AwLYlOsY2Wr07NnbF61cuvQTRo8aKXx/3/4DNt2U7YpLbyn3X4xklbx4ERInjEd6xi/w9/e39nCIBKxR87Ap25CbN2/CycnJ8HVxcTEi2rdDeXk5zn/3vRVHRnQba7TqePrChsQOHojGjQMRFhaGvLw8bN60CRcv/oC16zZYe2hEAFijlsCmbEO6d++BD9aswZbNm1BeXo6//KUVNmzchEGDB1t7aEQAWKOWwNMXREQS4eeUiYgkUm2nL5YvW4qFCxdAn5WF0NBQLHx3Edq373Df9SoqKnD16m9wcXHlPX3JJIqioLCwAL6+foYb55iCNUo1xawarY5bz23YuElxcHBQVry/SjmTdlYZPeZZxd3dXbn879/uu256RqYCgC++zH6lZ2SyRvmS+mVKjVbLOeVHO3dE23bt8e6ixQBuH1mENA3CCy++hMSkifdcNy8vD7oHGt1zGSI1+mvXodVqTVqWNUrWYEqNWvyccmlpKU6fPo2oqOg7G7G3R9eoaKSeOCEsX1JSgvz8fMOrsNC2H+VC1mPqqQTWKFmLKTVq8aack5OD8vJyeOu8jXKdtzf0+ixh+blzZsPby9PwahYcZOkhERlhjZLMrP7pi6SJk5Cdc8PwSs/ItPaQiIywRqkmWfzTF15eXqhXrx6y9dlGuT47Gzqdj7C8RqOBRqOx9DCI7oo1SjKz+JGyg4MDwsPDcfjwIUNWUVGBI4cPISIy0tKbIzIba5RkVi2fU45PeBXPjhmF8PC2aN++PZKTF6OoqAgjRj5THZsjMhtrlGRVLU150ODByMm5hhkzpkGflYWwsDDs+nSPTd/jlGoX1ijJSrp7X+Tn58Pby/P+CxL9SXbODbi5uVX7dlijVFmm1KjVP31BRER3sCkTEUmETZmISCJsykREEmFTJiKSCJsyEZFE2JSJiCTCpkxEJBE2ZSIiibApExFJhE2ZiEgibMpERBJhUyYikgibMhGRRNiUiYgkUi03uSei2s/Ly0vI5s6dq7rsX//6VyELCQkRsk8++UTIFi1aJGRffvmlKUO0STxSJiKSCJsyEZFE2JSJiCTCpkxEJBE+OLUaxcXFCVl0dLSQubi4CFmrVq2EzMfHR3U7V65cEbIvvvhCyE6dOiVkO3fuFLJLly6pbkd2fHBq9enVq5eQbdiwQcg8PDyqtB07OzshKyoqErL58+cL2fTp06u07ZrAB6cSEdkYNmUiIomwKRMRSYRNmYhIIpzos5C9e/cKWY8ePYTs8uXLQnbx4kUhy8jIEDK9Xq+67Y4dOwpZo0aNhKxt27ZCVlZWJmSTJk0SsnfeeUd12zLhRJ/5AgIChGz8+PFCNnLkSCFzd3cXsm3btqlup6CgQMjS09OFrG/fvkLWoUMHIVOb3O7cubOQ/frrr6rjsRZO9BER2Rg2ZSIiibApExFJhE2ZiEgivHWnhTRv3lzI7O3Ff/PUJkISExMtPp4GDRoIWdOmTYVs4MCBQqZ2ZdS1a9eEbOPGjZUcHVmD2gTetGnThKxhw4Ymvd+NGzeEbOjQoarLqk0oq1mwYIGQ7du3T8jUbgU6fPhwIZs1a5ZJ25UJj5SJiCTCpkxEJBE2ZSIiibApExFJhFf0Wcjy5cuF7Pnnnxeyxx57TMjUbrNpTampqULm5OQkZKGhoTUxHJPxir47mjRpImSZmZlCduvWLSFTe/6d2u1cBw8eLGRqt/i823uaqmXLlkL23XffCVlKSoqQDRgwoNLbrQ68oo+IyMawKRMRSYRNmYhIImzKREQSMfuKvmPHPsc7CxbgzJnTuHr1Kj7eug39+vc3fF9RFMyYPg0frFmN3NxcdOzUCcnJ7yFE5Yq32kTtmWEjRowQMkdHx5oYTpV89tlnQhYbG2uFkVQOaxR4+OGHhUxRxDn9NWvWCNnYsWNN2sbkyZOF7ObNmyatW1Vq+6J2lZ8tMvtI+feiIjwUGop3FyWrfn/B/HlY+t4SJC9ZimNffAXnhs7o27c3iouLqzxYIlOwRsmWmX2k3PPxXuj5uPrHXhRFwZLkxZj02uuI6dcPALD6g7UIDPDDrpQUDFY52iopKUFJSYnh64KCfHOHRGSENUq2zKLnlDMyMpCVlYWoqGhDptVq0b5DB6SmnlBdZ+6c2fD28jS8mgUHWXJIREZYoyQ7izZlvT4LAOCt0xnlOm8d9FlZquskTZyE7Jwbhld6RqYlh0RkhDVKsrP6rTs1Gg00Go21h1Flalc8paWlCVnPnj2F7MCBA9UxpEpr3bq1tYcgldpSo2rs7Owqve7dnhlZFc7OzkK2YcMGk9bNuss/qrbGokfKOp0PACD7T78sfbYeOh8fS26KqFJYoyQ7izbl4OBg+Pj44PDhQ4YsPz8f33z9NSIiIi25KaJKYY2S7Mw+fVFYWIj0//mvemZmBr5NS4OHpycCAwPxclw8Zs+aiZCQ5ggKCsL0aVPh6+dn9DlRourEGiVbZnZTPnXqJHp272b4OilxAgBg+NMjsGr1GoyfkIiioiK89OL/ITc3F506d8bu3Xts4qIJqh1Yo2TLeOvOajRlyhQhGzVqlJCpTf6pTRxWBx+V86jnzp0Tsi1btghZXFxctYypsnjrzjsiIiKEbMeOHULm5eUlZG+//baQqT3rrrS0tJKju83X11fI1Cb1unbtatL7rVy5UsheeOEF8wdWjXjrTiIiG8OmTEQkETZlIiKJsCkTEUmEE33VqEOHDkKm9jy+FStWCFlNTaK1a9dOyI4fPy5k0dHRQvb5559Xy5gqixN99xYSEiJkBw8eFLLGjRsLmdrVcmrPm7zbBLXae6ptu7nK7VPVbtO5e/duIRs0aJCQVXUy0tI40UdEZGPYlImIJMKmTEQkETZlIiKJcKKvhs2dO1fI2rRpI2QxMTFCVl5eXqVtq129t2nTJiHz8PAQskceeaRK264JnOgzn7+/v5DNmzdPyAYOHChk169fF7K1a9eqbqe/yn1FWrRoIWRqtxLduXOnkA0bNkzIaur5gFXBiT4iIhvDpkxEJBE2ZSIiibApExFJhE2ZiEgiVn9wal2zfv16IVO73HTq1KlC9uabb1Zp22PGjBGyLl26mJRR7XTlyhUhGzp0qJB99dVXQrZ48WIhS0pKqtJ45syZI2QzZswQMlv4pEVl8UiZiEgibMpERBJhUyYikgibMhGRRDjRV8POnz8vZG+99ZaQqd1PedmyZarvefXqVSELCwsTMrWJQrVLWGW7TzJZX3p6usXfc+bMmUL2xhtvWHw7toZHykREEmFTJiKSCJsyEZFE2JSJiCTCiT4JLFmyRMi8vb2FTO1ezAAwfvx4Idu4caOQpaSkCJna1VtEfzZu3DghU7v3sTn8/PyqtH5txSNlIiKJsCkTEUmETZmISCJsykREEuGDUyU1atQoIbvbFX1lZWVCdvLkSSHr16+fkBUUFFRidHLig1MtY8CAAUK2ZcsWIbO3F4/p1G4PCwCJiYlC1qpVKyGrV6+eKUO0WXxwKhGRjWFTJiKSCJsyEZFE2JSJiCTCK/ospGHDhkIWEREhZH369DHp/Vq3bi1kDg4OqsumpaUJWUxMjJAVFhaatG2q29Tqtn59sVWoPTtP7RmUgPqEstrk4aBBg4Rs69atqu9ZW/FImYhIImzKREQSYVMmIpKIWU157pzZ6NwxEl6e7mjs74tBTz6BHy9eNFqmuLgYCfFx8PPxRiMPLYYMHgS9Xm/RQRPdDWuUbJ1ZV/TF9O2NQYNj0a5tO9y6dQtvvjkFFy5cQNq35+Ds7AwAiHv5Jfzzn3vx/qrV0Gq1eDUhHnb29jhy9JhJ27DVq6U+/PBDIXvqqaeELCcnR8iysrKE7Oeffxaymzdvqm57yJAhQubk5CRkxcXFquvXFtk5NzBs6BDWqBlcXV2FTG3iOCgoSMg6deokZKmpqSZvOzs726T11SatbZUpV/SZ9emL3Z/uNfr6/VVr0NjfF6dPn8Kjjz6GvLw8rP1gDdat34iuXaMAACvfX42w0DZITT2BiIhIM3eByDysUbJ1VTqnnJ+XBwDw9Lh91HD69CmUlZUhKjrasEyLli3RODAQqSdOqL5HSUkJ8vPzDa+CgvyqDInICGuUbE2lm3JFRQUmTBiHjp06oXWbNgAAfZYeDg4OcHd3N1pW5+0Nvcp/0YHb5wC9vTwNr2bBQZUdEpER1ijZoko35YT4OFy4cAEbNm6q0gCSJk5Cds4Nwys9I7NK70f0B9Yo2aJKXdH3SkI89u7dgwMHDyMgIMCQ63x0KC0tRW5urtGRiD47GzofH9X30mg00Gg0lRmGVFxcXIRMbfJP7RaGahN9arRarWoeGxsrZL169RKyHTt2mLSd2oA1ahp/f38hCw4ONmldcyb11Kg94y8ykuf0zTpSVhQFryTEY1fKTuzbt1/45YWHt0WDBg1w+NAhQ/bjxYv49fJlRPCHTTWANUq2zqwj5YT4OHy0ZTO2btsOF1dXwxGeVquFk5MTtFotnhk1GklJE+Dh6QE3NzeMeyUBkZGRnNWmGsEaJVtnVlNeuWI5AKBHt2jjfNVqjBgxEgAwb/4C2Nvb46nYwSgpKUH37j2wKHmJhYZLdG+sUbJ1ZjXl4tJb913G0dERixYnY9Hi5EoPiqiyWKNk63jrTgvZvHmzkK1cuVLItm/fLmSmTsDl/fczt3+WkpIiZAMHDqz0dqjuKC0tFbKSkhIhU7ttbJcuXYRMbcIbADp37ixkale2nTp1SnX9uoQ3JCIikgibMhGRRNiUiYgkwqZMRCQRTvRZiNrzxsLDw4Vs9uzZQnbmzBkhy8zMNHnb69atE7KxY8eavD7VXWq3iD158qSQqd2m89D/XIDzB0Ux+U7Aql588cUqrV8b8EiZiEgibMpERBJhUyYikgibMhGRRDjRV42mTZsmZNevXxeyvXv3CpnahMfdbpUYFRUlZEePHjVhhESiVatWCVm7du2EzJzbmapNCk6ZMkXI1J4PWNfwSJmISCJsykREEmFTJiKSCJsyEZFE7IpLb1XtEhwLy8/Ph7eXp7WHUaNatGghZGvXrhWyiIgI1fWPHz8uZMOGDRMyc64StEXZOTdUbwdpaXWxRskyTKlRHikTEUmETZmISCJsykREEmFTJiKSCK/ok8DFixeFrGPHjlYYCRFZG4+UiYgkwqZMRCQRNmUiIomwKRMRSYRNmYhIImzKREQSYVMmIpIImzIRkUSka8qKItVN68iG1FTtsEapskypHemacmFhgbWHQDaqpmqHNUqVZUrtSHc/5YqKCly9+htcXFxRWFiAZsFBSM/IhKtr9d8nt7oUFOTXiv0A5NwXRVFQWFgAX18/2NtX/3EGa1RuMu6LOTUq3b0v7O3t4e8fAACws7MDALi6utXIzcurW23ZD0C+fdFqtTW2LdaobZBtX0ytUelOXxAR1WVsykREEpG6KWs0Gkye8gY0Go21h1IltWU/gNq1L5ZQW34etWU/ANvfF+km+oiI6jKpj5SJiOoaNmUiIomwKRMRSYRNmYhIImzKREQSkbYpL1+2FA82bwatqzMe7dwR33zztbWHdF/Hjn2OJ/7eH8FNGsPRoT52paQYfV9RFEyfNhVBgQFwd3NBr8d74NJPP1lptHc3d85sdO4YCS9PdzT298WgJ5/Aj3964nZxcTES4uPg5+ONRh5aDBk8CHq93kojtg7WqPXU5hqVsilv/fhjJCVOwOQpb+BE6jd4KDQMMX16Izs729pDu6ffi4rwUGgo3l2UrPr9BfPnYel7S5C8ZCmOffEVnBs6o2/f3iguLq7hkd7bsWOfY+wLL+DzY19iz95/oexWGfr06YWioiLDMokTxmPPnk/x4eYt2H/wEK5e/Q2xgwdacdQ1izVqXbW5RqX8nPKjnTuibbv2eHfRYgC3bwAT0jQIL7z4EhKTJlp5dKZxdKiPj7duQ7/+/QHcPgIJbtIYCa+8ilfHjQcA5OXlITDAD++vWoPBsbHWHO49Xbt2DY39fbH/4CE8+uhjyMvLQ4CfD9at34gnnnwSAHDxhx8QFtoGR499gYiISCuPuPqxRuVSm2pUuiPl0tJSnD59GlFR0YbM3t4eXaOikXrihBVHVjUZGRnIysoy2i+tVov2HTogNVXu/crPywMAeHp4AgBOnz6FsrIyREXf2ZcWLVuicWCgTf+OTMUalU9tqlHpmnJOTg7Ky8vhrfM2ynXe3tDrs6w0qqr7Y+zeOp1RrvPWQZ8l735VVFRgwoRx6NipE1q3aQMA0Gfp4eDgAHd3d6Nldd7eUu+LpbBG5VLbalS6W3eSXBLi43DhwgUcOnzU2kMhUlXbalS6I2UvLy/Uq1cP2XrjCRN9djZ0Oh8rjarq/hh79p9mf/XZeuh85NyvVxLisXfvHuz77AACAgIMuc5Hh9LSUuTm5hotr8/OlnZfLIk1Ko/aWKPSNWUHBweEh4fj8OFDhqyiogJHDh9CRKS8J+fvJzg4GD4+Pkb7lZ+fj2++/lq6SQdFUfBKQjx2pezEvn37ERwcbPT98PC2aNCgAQ4furMvP168iF8vX7bp35GpWKPWV5trVMrTF/EJr+LZMaMQHt4W7du3R3LyYhQVFWHEyGesPbR7KiwsRPqlS4avMzMz8G1aGjw8PREYGIiX4+Ixe9ZMhIQ0R1BQEKZPmwpfPz/D7LcsEuLj8NGWzdi6bTtcXF2R9d9zcFqtFk5OTtBqtXhm1GgkJU2Ah6cH3NzcMO6VBERGRkr3l7e6sEatqzbXqJQfiQOAZUvfw8KFC6DPykJYWBgWvPMuOnSIsPaw7uno0SPo2b2bkA9/egRWrV4DRVEwY/o0rFm9Crm5uejUuTMWL16C5g8+aIXR3p2jg/q/1StXrcaIESMB3P5g/sSkRHz80RaUlJSge/ceWJS8BD6S/9fQklij1lOba1TapkxEVBdJd06ZiKguY1MmIpIImzIRkUTYlImIJMKmTEQkETZlIiKJsCkTEUmETZmISCJsykREEmFTJiKSCJsyEZFE/h8H62qR2GW6/wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0439\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWUAAADOCAYAAADmKRMiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAaYElEQVR4nO3deVhU9f4H8DdoDCowbDLsimuLQUEI6q9uopgLmpVLi0suXTMDXHDJLJfMXVMxl1RSM7XFfbmZCyrVVbsiaty0NAw3BtEYFgMUzu8Pa2zu96QDDMx38P16nnmeeHNmzueMHz+e5jtzxq6o5JYCIiKSgr21CyAiojs4lImIJMKhTEQkEQ5lIiKJcCgTEUmEQ5mISCIcykREEuFQJiKSCIcyEZFEOJSJiCTCoVwOa9ashqNDbRw79h9rl1KlPlq2FC+/2BtNGgfB0aE2Bg8aaO2SyEz3Q49euHABU9+bgv9rHQlvL0/4+egQ3T4K+/bttXZpFsGhTII5c2bjwIFkPPzww6hdu7a1yyEysWP7NsydMxuNGjfBpMlT8Nb4t1GQX4AunTpi9epV1i6v0vg3jgR79u5HYGAg7Ozs4OGmtXY5RCb+8Y+n8fO5DHh6ehqz1/45BC2fCMOUyZPQv/+r1ivOAnimXEmDBw2Eh5sWmZmZeK57N3i4adGoYSCWLlkMAPjh1Ck806E93F1d0LRJI2xYv97k/tevX8e4saMR9vhj8HDTor6HG7p17YKTJ04I+/r111/xwnPd4e7qggA/H4xOGIk9X++Go0NtHDx4wGTbo0ePoGtMZ3h5usNN64z27driu+++NeuYGjRoADs7u4o9ISSdmtajDz/yiMlABgCNRoOOnTri0sWLyM/PL+czJBcOZQsoLS3Fs11j4O8fgGnTZyCwQQMMj4/DmjWr0a1rF4SFhuH9adPh7OSMQQNfRUZGhvG+GRm/YNu2bejUuTNmzp6DESNHIf2HHxDdPgqXL182bldYWIiOz0Rj//59eGPYmxg77i0c/ve/8fb48UI9ycn70T6qLfLz8vD2hHcw5b2pMOQa0LFDNL7//mi1PCckl/uhR7Oy9Khbty7q1q1boftLo6jklsKbebePVqxUACjf/vuwMevTt58CQJny3lRjlpWdo9SpU0exs7NTPlm7zpifPJWuAFDenvCOMTPkFyo3ikpM9nP6p7OKRqNR3p04yZjNnDVbAaB88eUmY5abV6A0b/6gAkDZvWevUlRyS/m9+KbSpElTJbpDB+X34pvGbX8z5CsNg4KUdu3bl+uY69Wrp/Tp28/qzz1v7NG73dL/e1pxdHRUXn6lj9X/DCp745myhQwYOMj4366urmjWrDnq1auHHj17GvNmzZvD1dXV5CxEo9HA3v72H0NpaSmuXbsGJycnNGvWHGnHjxu3+/rr3fDz80NM167GzNHREQMH3dkvAJxIS8PZsz+jd++XcO3aNeTk5CAnJweFhYVo2zYK36SkoKyszOLHT/KrqT1648YNvPzSi6hTpw6mvj/N/CdEUlzoswBHR0fUr1/fJHPRusDPz194bdZFq0Xub78Zfy4rK8OixIVYtmwpzmdkoLS01Pg7dw93439n/pqJoEaNhMdr3LiJyc9nz54FAAweNOBv6zUYDHBzczPz6KgmqKk9Wlpair6vvIwff/wvtm7fAV9f33veR3YcyhZQq1atcuUK7nwD18wZ0zF50kT0f3UAJk6cDDd3N9jb22P0qFEVOqP98z7TZ8xEcEiI6jZOTk7lflyybTW1R4e+PgS7du3EqtWfoG3bqHLXIiMOZSvbvGkT/vH001j20XKT3GDIhYenh/HnwAaBOP3jj1AUxeRM5Ny5syb3a9S4EQDA2cUF7dq1r8LK6X4ha4++NW4M1qxehTlz56H3iy9W+HFkw9eUraxWrVpQFNPvrt345Ze4dOmSSRYd3QGXLl3Cju3bjVlRURGSVq402S40NAyNGjfG/HnzUFBQIOzv6tWrFqye7gcy9ui8uXPwwbx5GDN2HN6MjSvP4UiPZ8pW1qlzZ0x7fypeGzwIka1aIf2HH7Bh/ToENWpkst3g1/6JJUsWo1/fVzDszVj4+Phgw/p1cHR0BADjmYm9vT2WLF2GZ7vG4PHHgtGvX3/4+vnh8qVLOHjwAFycXbBpy9a71rRzx3acPHkSAHDz5k38cOoUpk97HwAQE9MVjwYHW/ppIInJ1qNbt2zB+LfGoUmTpnjwwYew7tNPTX7frn176HQ6Cz8L1YdD2crGjnsLN24U4rMNG/DlF5/jsccfx+at2zDhbdP3djo5OeGr3XswckQ8PlyUCCcnJ7zSpw8iI1vjxd49jY0P3P7E08FD32D6tPexdMliFBQUQOftjZbhLTHotdfuWdPmzZux9pM1xp/T0o4jLe32Krufvz+H8n1Gth49efL2h1bOnv0ZAwf0F36/e89emx7KdkUlt5R7b0aySly4AKMTRuFcxq/w8/OzdjlEAvZo+XAo25Dff/8dderUMf5cVFSEiPAnUFpaih/++6MVKyO6jT1aeXz5wob07tUDAQGBCAkJgcFgwPp163DmzGmsWv2JtUsjAsAetQQOZRsSHd0BHyclYcP6dSgtLcVDDz2MT9auQ89evaxdGhEA9qgl8OULIiKJ8H3KREQSqbKXL5YuWYx58+ZCn5WF4OBgzJu/AOHhLe95v7KyMly5chlOTs68pi+ZRVEUFBTkw8fH13jhHHOwR6m6lKtHq+LSc5+sXac4ODgoy5avUI6nnVQGDhqsuLq6KpkXL9/zvucyzisAeOOt3LdzGefZo7xJfTOnR6vkNeUn27RC2BPhmL9gIYDbZxZNGjXE0DeGYfSYsXe9r8FggK6+x123IVKjv3oNWq3WrG3Zo2QN5vSoxV9TLikpQWpqKqKi2t3Zib092ka1w5HDh4Xti4uLkZeXZ7wVFNj2V7mQ9Zj7UgJ7lKzFnB61+FDOyclBaWkpvHReJrnOywt6fZaw/ayZM+Dl6W68NQ5qaOmSiEywR0lmVn/3xZix45Cdc914O5dx3tolEZlgj1J1svi7Lzw9PVGrVi1k67NNcn12NnQ6b2F7jUYDjUZj6TKI/hZ7lGRm8TNlBwcHhIaGIjl5vzErKyvDgeT9iIiMtPTuiMqNPUoyq5L3KcfFj8DgQQMQGhqG8PBwJCYuRGFhIfr1f7UqdkdUbuxRklWVDOWevXohJ+cqpkyZBH1WFkJCQrBtx06bvsYp1SzsUZKVdNe+yMvLg5en+703JPof2TnX4eLiUuX7YY9SRZnTo1Z/9wUREd3BoUxEJBEOZSIiiXAoExFJhEOZiEgiHMpERBLhUCYikgiHMhGRRDiUiYgkwqFMRCSRKvviVCKqORwcHIQsLi5OyGJiYlTvv2XLFiGbP39+ZcuqkXimTEQkEQ5lIiKJcCgTEUmEQ5mISCJc6COiexo1apSQTZ061ez7Dxw40JLl1Gg8UyYikgiHMhGRRDiUiYgkwqFMRCQRLvQRkQm1T9r17dvXrPuuWrVKNf/ll18qUdH9hWfKREQS4VAmIpIIhzIRkUQ4lImIJGJXVHJLsXYRf5WXlwcvT3drl2ERISEhQhYbGytkPXr0EDIXFxchs7OzU92PwWAQsi+//FLIEhMThezEiROqj2mLsnOuqz5vllaTetTR0VHIMjMzhczDw6NS+7lx44ZZ261du1bIpkyZImRXrlypVD3WYk6P8kyZiEgiHMpERBLhUCYikgiHMhGRRLjQVwEajUbIunfvLmRJSUlCprawoubixYtCFhAQoLqtopj3R6i22PLuu+8K2QcffGDW48mGC33lt3r1aiF7+eWXhWz58uVC9thjj5mVAep9b27fFhcXC9krr7wiZJs3bzbr8ayJC31ERDaGQ5mISCIcykREEuFQJiKSCBf6KiAsLEzIjh49atZ9k5OThWzWrFlCduTIESF76KGHzNoHAMTHxwtZr169hKywsFDIXnrpJSHbuXOn2fu2Fi703V2LFi2E7LvvvhOy9957T8hmz55t1j7+bqGvQ4cOZmVt27Y1az9qfduqVSshS09PN+vxqgsX+oiIbAyHMhGRRDiUiYgkwqFMRCSRci/0paQcwgdz5+L48VRcuXIFn3+xEd2efdb4e0VRMGXyJHyctBK5ublo1bo1EhM/RJOmTc16fFtYRPnss8+ErFu3bkKm9imoESNGCFlpaallCvuL2rXFr19UuwTi2LFjhSw1NVXIwsPDLVNYFfpzEYU9qn6Z148//ljI1D695+fnJ2RXr161TGF/4eDgIGRDhgwRMrXvDFSjtjjeqVMnIVO71G11qZKFvhuFhXg0OBjzF4jX5gWAuXNmY/GHi5C4aDFSvvkO9erWQ0xMZxQVFZV3V0QVwh4lW1bub7N+pmMnPNNR/NcHuH0GsihxIca9NR5d/zhzXPnxKgT6+2Lb1q3o1bu3cJ/i4mKTz7bn5+eVtyQiE+xRsmUWfU05IyMDWVlZiIpqZ8y0Wi3CW7bEkSOHVe8za+YMeHm6G2+NgxpasiQiE+xRkp1Fh7JenwUA8NLpTHKdlw76rCzV+4wZOw7ZOdeNt3MZ5y1ZEpEJ9ijJrtwvX1iaRqNRvRSmzKKiooQsPz9fyOLi4qqjHFW3bt0SMrVPb6lp0KCBpcuxabbYo15eXkLWt29fIduwYYOQVcWinpqSkhIhU/seSTVqi38RERFClpCQIGTvvPOOWfuwFoueKet03gCAbL3eJNdn66Hz9rbkrogqhD1KsrPoUA4KCoK3tzeSk/cbs7y8PHx/9CgiIiItuSuiCmGPkuzK/fJFQUEBzp09a/z5/PkMnEhLg5u7OwIDA/FmbBxmTJ+GJk2aomHDhpg8aSJ8fH1N3idKVJXYo2TLyj2Ujx37D56Jbm/8eczo26/Z9OnbDytWJmFUwmgUFhZi2BuvIzc3F63btMH27TvN/hokospij5It46U778Hf31/I1C4HqPY9YmqLLZXh5OSkmpeVlQlZcHCwkK1du1bIgoKChOzatWtCZuljqQq8dOcdQ4cOFbJFixYJmdrlXDdu3FglNVWUq6urkJ06dUrIfH19hUxtIVPt+/2qCy/dSURkYziUiYgkwqFMRCQRDmUiIolY/RN9srt48aKQqX0SSe0yhD4+PkJ25cqVCteyfft21Vz3Px8ZBtQX5tzc3MzaT1JSUvkKI+nExMQIWVpampDt2LGjGqqpnNzcXCFLSUkRst4qF5Pq2rVrVZRUpXimTEQkEQ5lIiKJcCgTEUmEQ5mISCJc6KuATz/9VMhiY2OFbMKECWZtp/aJPDV5eerfePHUU0+ZdX81+/btE7Lx48dX+PGo+kVGihdSio6OFrJjx44JmdonUW2B2ncQqmVql9SVHc+UiYgkwqFMRCQRDmUiIolwKBMRSYRDmYhIInz3RQWoXaN1yJAhQvb6668Lmdo1kWfOnClkv/zyi5ANGzZMtR6DwSBk5l4zVm3f5r4bhOTQokULIatVq5aQKYpUl06vFLVjUcts8ZIBPFMmIpIIhzIRkUQ4lImIJMKhTEQkES70VcDhw4eFbPjw4UK2ePFiIevTp49Z2fHjx4VM7frMAODt7a2am0NtP0RkPTxTJiKSCIcyEZFEOJSJiCTCoUxEJBEu9FnIsmXLhGz37t1Ctn79eiGLiIgQstDQUCFTu14sYP4ntdS+OLOwsNCs+5LtCwgIMCu7cOFCdZRjtgYNGghZjx49hOzkyZNCNmfOnCqpqSrxTJmISCIcykREEuFQJiKSCIcyEZFEuNBXhc6fPy9kTz/9tJDVr1+/UvtJSEgQMrUvaFVbwLHVL86kO7Kzs4VM7fKrvr6+Qqa2oGzNhT4vLy8h27Jli5CpXZp03rx5QqZ2WVvZ8UyZiEgiHMpERBLhUCYikgiHMhGRRLjQV83UFtYuXrxYqcfMz8+v1P3Jtm3btk3IioqKhKxu3bpCtnz5ciE7dOiQkP32228VrK581Batg4ODhWzTpk1Ctnnz5iqpqbrxTJmISCIcykREEuFQJiKSSLmG8qyZM9CmVSQ83V0R4OeDni88j5/OnDHZpqioCPFxsfD19oKHmxYv9uoJvV5v0aKJ/g57lGxduRb6UlIOYcjQoXgi7AncunUL7747AV26dELaiVOoV68eAGB0wij861+78On6DdBqtRgRH4fevXrgwMGUKjkAor9ij96m9v2QaotoHh4eQjZy5Egh27Fjh5CdOHFCyNQWGAH1TxNu3LhRyMLCwoTs8uXLQta3b1+z921ryjWUt+/YZfLz8hVJCPDzQWrqMTz55FMwGAxY9XESVq9Zi7ZtowAAHy1fiZDgFjhy5DAiIiItVzmRCvYo2bpKvaac98fnyt3d3AEAqanHcPPmTUS1a2fcpvmDDyIgMBBHVL4BGrj9FrG8vDzjLT8/rzIlEZlgj5KtqfBQLisrQ0LCSLRq3RqPtGgBANBn6eHg4ABXV1eTbXVeXtBnZak+zqyZM+Dl6W68NQ5qWNGSiEywR8kWVXgox8fFIj09HZ+sXVepAsaMHYfsnOvG27mM85V6PKI/sUfJFlXoE33D4+Owa9dO7N2XDH9/f2Ou89ahpKQEubm5Jmci+uxs6Ly9VR9Lo9FAo9FUpAyiv3W/96jaZSzVFvrUjB8/3qwsMzNTyPbu3av6mIMHDxYytcuLqn1y8LnnnhOymrKop6ZcZ8qKomB4fBy2bd2C3bv3ICgoyOT3oaFheOCBB5C8f78x++nMGVzIzEREJBdQqOqxR8nWletMOT4uFp9tWI8vNm6Ck7Mzsv54DU6r1aJOnTrQarV4dcBAjBmTADd3N7i4uGDk8HhERkZyVZuqBXuUbF25hvJHy5YCADq0b2ear1iJfv36AwBmz5kLe3t7vNS7F4qLixEd3QELEhdZqFyiu2OPkq0r11AuKrl1z20cHR2xYGEiFixMrHBRRBXFHiVbx0t3EtVAapeIVftknNon7cwVGBgoZAMHDlTd9tSpU0K2YMECIfv222+F7PTp0xWoznbxgkRERBLhUCYikgiHMhGRRDiUiYgkwoW++4jaoo6Dg4OQlZSUVEc5VIVyc3OFrEuXLkL21VdfCZlOpxOyffv2CVl6errZ9UyaNEnIDH9cLIpM8UyZiEgiHMpERBLhUCYikgiHMhGRRLjQdx8JDQ0VMicnJyG7fv16dZRD1ezkyZNCVplP9FHV4JkyEZFEOJSJiCTCoUxEJBEOZSIiiXChrwa4ceNGhe+r9v1nK1eurEw5RFQJPFMmIpIIhzIRkUQ4lImIJMKhTEQkES701QBq33UWEREhZI8//riQHTx4sEpqIqKK4ZkyEZFEOJSJiCTCoUxEJBEOZSIiiXChrwZQ+0Rf9+7dq78QIqo0nikTEUmEQ5mISCIcykREEpFuKCuKYu0SyEZVV++wR6mizOkd6YZyQUG+tUsgG1VdvcMepYoyp3fsikpuSfXPfllZGa5cuQwnJ2cUFOSjcVBDnMs4D2dnF2uXVmH5+Xk14jgAOY9FURQUFOTDx8cX9vZVf57BHpWbjMdSnh6V7i1x9vb28PPzBwDY2dkBAJydXeDiIseTWxk15TgA+Y5Fq9VW277Yo7ZBtmMxt0ele/mCiOh+xqFMRCQRqYeyRqPB2xPegUajsXYplVJTjgOoWcdiCTXl+agpxwHY/rFIt9BHRHQ/k/pMmYjofsOhTEQkEQ5lIiKJcCgTEUmEQ5mISCLSDuWlSxajWdPG0DrXw5NtWuH7749au6R7Skk5hOe7P4ugBgFwdKiNbVu3mvxeURRMnjQRDQP94erihE4dO+Dszz9bqdq/N2vmDLRpFQlPd1cE+Pmg5wvP46czZ0y2KSoqQnxcLHy9veDhpsWLvXpCr9dbqWLrYI9aT03uUSmH8heff44xoxPw9oR3cPjI93g0OARdu3RGdna2tUu7qxuFhXg0OBjzFySq/n7unNlY/OEiJC5ajJRvvkO9uvUQE9MZRUVF1Vzp3aWkHMKQoUNxKOVb7Nz1FW7euokuXTqhsLDQuM3ohFHYuXMHPl2/AXv27ceVK5fRu1cPK1Zdvdij1lWTe1TK9yk/2aYVwp4Ix/wFCwHcvgBMk0YNMfSNYRg9ZqyVqzOPo0NtfP7FRnR79lkAt89AghoEIH74CIwYOQoAYDAYEOjvi+UrktCrd29rlntXV69eRYCfD/bs248nn3wKBoMB/r7eWL1mLZ5/4QUAwJnTpxES3AIHU75BRESklSuueuxRudSkHpXuTLmkpASpqamIimpnzOzt7dE2qh2OHD5sxcoqJyMjA1lZWSbHpdVqEd6yJY4ckfu48gwGAIC7mzsAIDX1GG7evImodneOpfmDDyIgMNCm/4zMxR6VT03qUemGck5ODkpLS+Gl8zLJdV5e0OuzrFRV5f1Zu5dOZ5LrvHTQZ8l7XGVlZUhIGIlWrVvjkRYtAAD6LD0cHBzg6upqsq3Oy0vqY7EU9qhcalqPSnfpTpJLfFws0tPTsT/5oLVLIVJV03pUujNlT09P1KpVC9l60wUTfXY2dDpvK1VVeX/Wnv0/q7/6bD103nIe1/D4OOzatRO7v94Lf39/Y67z1qGkpAS5ubkm2+uzs6U9Fktij8qjJvaodEPZwcEBoaGhSE7eb8zKyspwIHk/IiLlfXH+XoKCguDt7W1yXHl5efj+6FHpFh0URcHw+Dhs27oFu3fvQVBQkMnvQ0PD8MADDyB5/51j+enMGVzIzLTpPyNzsUetryb3qJQvX8TFj8DgQQMQGhqG8PBwJCYuRGFhIfr1f9Xapd1VQUEBzp09a/z5/PkMnEhLg5u7OwIDA/FmbBxmTJ+GJk2aomHDhpg8aSJ8fH2Nq9+yiI+LxWcb1uOLjZvg5OyMrD9eg9NqtahTpw60Wi1eHTAQY8YkwM3dDS4uLhg5PB6RkZHS/eWtKuxR66rJPSrlW+IAYMniDzFv3lzos7IQEhKCuR/MR8uWEdYu664OHjyAZ6LbC3mfvv2wYmUSFEXBlMmTkLRyBXJzc9G6TRssXLgITZs1s0K1f8/RQf3f6o9WrES/fv0B3H5j/tgxo/H5ZxtQXFyM6OgOWJC4CN6S/6+hJbFHracm96i0Q5mI6H4k3WvKRET3Mw5lIiKJcCgTEUmEQ5mISCIcykREEuFQJiKSCIcyEZFEOJSJiCTCoUxEJBEOZSIiiXAoExFJ5P8B+NsdNpMSQecAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9008\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWUAAADOCAYAAADmKRMiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAZc0lEQVR4nO3deVxU9foH8A9gDCgwbM6wC4qaZVogivrzliDu2m1xueW+lFaCKahl5XLLXW+KW6XmUmKL5YJ0zQSN6rokmek1FcMslUH0wjAooHB+f1Cj0znqDAzMd+Dzfr3O69U8c5bnjI+P38535hyHkrKbEoiISAiOtk6AiIhuYVMmIhIImzIRkUDYlImIBMKmTEQkEDZlIiKBsCkTEQmETZmISCBsykREAmFTJiISCJuyBTZu3AAX5wY4cuR7W6dSo959ZzWeGTwI4c3C4OLcAGNGj7J1SmSm+lCjv/32G97852z8X6do+Gl8EeivRVy3GOzd+5WtU7MKNmWSWbRoIfbty8ADDzyABg0a2DodIhOpO3dg8aKFaNosHDNnzcYrr06HociAPr16YsOG9bZOr9r4N45k9nyVjpCQEDg4OMDHS23rdIhMPProYzhzNge+vr7G2Njnnkf7dpGYPWsmhg8fYbvkrIAj5WoaM3oUfLzUOH/+PJ74e3/4eKnRNDQEq1etBAAc/+kn9OjeDd6eHmge3hRbUlJMtr969SqmTU1C5CMPw8dLjcY+Xujfrw+O/fij7Fi//vornnri7/D29EBwoD+SEidhz5e74eLcAPv37zNZ99Chg+jXtzc0vt7wUrujW2xXfPfdt2adU5MmTeDg4FC1D4SEU9dq9IEHHzRpyACgUqnQs1dPXPj9dxQVFVn4CYmFTdkKysvL8Xi/vggKCsacufMQ0qQJJibEY+PGDejfrw8iIyLx1py5cHdzx+hRI5CTk2PcNifnF+zYsQO9evfG/IWL8PKkyThx/DjiusXg4sWLxvWKi4vRs0cc0tP34oUXX8LUaa/gwH/+g+mvvirLJyMjHd1iuqJIr8f0117H7H++icKCQvTsHofDhw/VymdCYqkPNZqbq0PDhg3RsGHDKm0vjJKymxIX85Z316yVAEjf/ueAMTZk6DAJgDT7n28aY7l5+ZKrq6vk4OAgbfpgszF+7KcTEgBp+muvG2OFRcXStZIyk+P8fDpbUqlU0hszZhpj8xcslABIn3z6mTFWoDdILVveLwGQdu/5SiopuyldL70hhYc3l+K6d5eul94wrvu/wiIpNCxMiu3WzaJzbtSokTRk6DCbf/ZcWKN3W07892fJxcVFeubZITb/M6juwpGylYwcNdr4356enmjRoiUaNWqEpwcMMMZbtGwJT09Pk1GISqWCo2PlH0N5eTmuXLkCNzc3tGjREkd/+MG43pdf7kZgYCD69utnjLm4uGDU6FvHBYAfjx5FdvYZDBr0D1y5cgX5+fnIz89HcXExunaNwTeZmaioqLD6+ZP46mqNXrt2Dc/8YzBcXV3x5ltzzP9ABMWJPitwcXFB48aNTWIeag8EBgbJrs16qNUo+N//jK8rKiqwPHkZ3nlnNc7l5KC8vNz4nrePt/G/z/96HmFNm8r216xZuMnr7OxsAMCY0SPvmG9hYSG8vLzMPDuqC+pqjZaXl2Pos8/g5Mn/YvvOVAQEBNxzG9GxKVuBk5OTRXEJt57ANX/eXMyaOQPDR4zEjBmz4OXtBUdHRyRNnlylEe2f28ydNx9t2rZVXMfNzc3i/ZJ9q6s1On7c80hL24X1Gzaha9cYi3MREZuyjX3+2Wd49LHH8M6775nECwsL4OPrY3wd0iQEP588CUmSTEYiZ89mm2zXtFlTAIC7hwdiY7vVYOZUX4hao69Mm4KNG9Zj0eIlGDR4cJX3IxpeU7YxJycnSJLps2u3fvopLly4YBKLi+uOCxcuIHXnTmOspKQE69auNVkvIiISTZs1w9tLlsBgMMiOd/nyZStmT/WBiDW6ZPEi/GvJEkyZOg0vTYi35HSEx5GyjfXq3Rtz3noTY8eMRnTHjjhx/Di2pGxGWNOmJuuNGfscVq1aiWFDn8WLL02Av78/tqRshouLCwAYRyaOjo5YtfodPN6vLx55uA2GDRuOgMBAXLxwAfv374OHuwc+27b9rjntSt2JY8eOAQBu3LiB4z/9hLlz3gIA9O3bDw+1aWPtj4EEJlqNbt+2Da++Mg3h4c1x//2tsPnDD03ej+3WDVqt1sqfQu1hU7axqdNewbVrxfhoyxZ8+snHePiRR/D59h14bbrpdzvd3Nzw7917MOnlBKxYngw3Nzc8O2QIoqM7YfCgAcbCByp/8bT/628wd85bWL1qJQwGA7R+fmgf1R6jx469Z06ff/45Pti00fj66NEfcPRo5Sx7YFAQm3I9I1qNHjtW+aOV7OwzGDVyuOz93Xu+suum7FBSdlO692okquRlS5GUOBlnc35FYGCgrdMhkmGNWoZN2Y5cv34drq6uxtclJSXoENUO5eXlOP7fkzbMjKgSa7T6ePnCjgwa+DSCg0PQtm1bFBYWImXzZpw69TPWb9hk69SIALBGrYFN2Y7ExXXH++vWYUvKZpSXl6NVqwew6YPNGDBwoK1TIwLAGrUGXr4gIhIIv6dMRCSQGrt8sXrVSixZshi63Fy0adMGS95eiqio9vfcrqKiApcuXYSbmzvv6UtmkSQJBkMR/P0DjDfOMQdrlGqLRTVaE7ee2/TBZsnZ2Vl657010g9Hj0mjRo+RPD09pfO/X7zntmdzzkkAuHCxeDmbc441ykXoxZwarZFryl06d0Rkuyi8vXQZgMqRRXjTUIx/4UUkTZl6120LCwuhbexz13WIlOguX4FarTZrXdYo2YI5NWr1a8plZWXIyspCTEzsrYM4OqJrTCwOHjggW7+0tBR6vd64GAz2/SgXsh1zLyWwRslWzKlRqzfl/Px8lJeXQ6PVmMS1Gg10ulzZ+gvmz4PG19u4NAsLtXZKRCZYoyQym3/7YsrUacjLv2pczuacs3VKRCZYo1SbrP7tC19fXzg5OSFPl2cS1+XlQav1k62vUqmgUqmsnQbRHbFGSWRWHyk7OzsjIiICGRnpxlhFRQX2ZaSjQ3S0tQ9HZDHWKImsRr6nHJ/wMsaMHomIiEhERUUhOXkZiouLMWz4iJo4HJHFWKMkqhppygMGDkR+/mXMnj0TutxctG3bFjtSd9n1PU6pbmGNkqiEu/eFXq+Hxtf73isS/UVe/lV4eHjU+HFYo1RV5tSozb99QUREt7ApExEJhE2ZiEggbMpERAJhUyYiEgibMhGRQNiUiYgEwqZMRCQQNmUiIoGwKRMRCYRNmYhIIGzKREQCYVMmIhIImzIRkUDYlImIBMKmTEQkEDZlIiKBsCkTEQmETZmISCBsykREAqmRp1lTpQULFshi/fv3l8VatGghi61bt87q+Zw5c0YW++KLL2SxY8eOWf3YRGQejpSJiATCpkxEJBA2ZSIigbApExEJhBN9VuLv7y+LPf/887KY0mRbq1atZLGRI0fKYqmpqYrHdnFxkcW0Wq0sNmPGDLNyVDr2/v37FY9NthccHGzWeoMHD5bF/Pz8rJpLRESEYvyxxx6TxXbu3CmLKf39qI73339fFjt+/LhVj2FtHCkTEQmETZmISCBsykREAmFTJiISiENJ2U3J1kncTq/XQ+Prbes07krpl3q9evWSxYqKimQxpQmPsrIyq+R1L0oTeGvXrpXFjhw5IotFRUXVSE7WlJd/FR4eHjV+nJqoUXd3d1nswQcflMXGjh0ri40YMUIWk6Sq/7V2cHCw6v5qYp/m7u/333+XxSIjI2WxK1euVDkXS5hToxwpExEJhE2ZiEggbMpERAJhUyYiEgh/0XcP48ePl8USExNlMYPBIIt17NhRFqutST0l5k7W5ebm1nAm9FdPPvmkLKY0CWttp0+flsXMvXWr0rYBAQGK67q5uVmW2G2aN28uiz388MNmbRsUFCSLubq6VjmX2sCRMhGRQNiUiYgEwqZMRCQQNmUiIoFYPNGXmfk1/rV4MX74IQuXLl3Cx59sRf/HHze+L0kSZs+aiffXrUVBQQE6duqE5OQVCFe4WG8PnnnmGbPW+/TTT2WxEydOWDsds3l6espisbGxZm0r+q0N70X0Gn3jjTdksaSkJKseY9euXbKY0sTc0qVLZTGlX8HZ0qOPPiqLpaenm7Wt0udw6dKlaudUkyweKV8rLsZDbdrg7aXJiu8vXrQQK1csR/Lylcj85js0atgIffv2RklJSbWTJTIHa5TsmcUj5R49e6FHT/l9HoDKEcjy5GWY9sqr6PfHU5vXvr8eIUEB2LF9OwYOGiTbprS0FKWlpcbXRUV6S1MiMsEaJXtm1WvKOTk5yM3NRUzMrf9NVqvViGrfHgcPHlDcZsH8edD4ehuXZmGh1kyJyARrlERn1aas01X+6EDzl0cRaTVa6O7wg4QpU6chL/+qcTmbc86aKRGZYI2S6Gz+iz6VSgWVSmXrNAAAY8aMkcU6d+4siyndNjAvL69GcqoqpWevKf0ySulclGL1mbVr1MnJSRZr2LChWdsqTVzNmjVLFlO6/aq92rRpkyxmbo1evnxZFisvL692TjXJqiNlrbayEeTpdCZxXZ4OWis/oJGoKlijJDqrNuWwsDD4+fkhI+PW11X0ej0OHzqEDh2irXkooiphjZLoLL58YTAYcDY72/j63Lkc/Hj0KLy8vRESEoKXJsRj3tw5CA9vjtDQUMyaOQP+AQEm3xMlqkmsUbJnFjflI0e+R4+4bsbXU5Iq75g2ZOgwrFm7DpMTk1BcXIwXXxiHgoICdOrcGTt37oKLi4v1sia6C9Yo2TM+o+82aWlpsliPHj1ksV9++UUWU3pGX/Zto7XaNn36dFls9uzZspjShMnQoUNlsQ8//NA6idUge3lGn9KkXu/evc3aVmmi7/r161XORTQ+Pj6ymNIvTDUajSx26tQpWaxLly6yWG09j08Jn9FHRGRn2JSJiATCpkxEJBA2ZSIigdj8F30iadeunVnrzZkzRxaz5aSeEnd3d7PWO3nypCy2fft2a6dDt7l27ZospnTr1/po3Lhxsljjxo3N2lZpAs+Wk3pVxZEyEZFA2JSJiATCpkxEJBA2ZSIigXCi7x4KCgpkscOHD9d+IhaKiooya72rV6/KYgaDwdrpEMlERkbKYkq/RDVXVlZWddIRBkfKREQCYVMmIhIImzIRkUDYlImIBMKmTEQkEH774jZ9+/aVxZS+nSDaT6qVPPTQQ7ZOgeiu4uLiZDFzH1Cr1+tlsRUrVlQ7JxFwpExEJBA2ZSIigbApExEJhE2ZiEggnOi7zaFDh2ydQpU899xzspivr69Z265fv97K2RDJtW7dWhYbP368LCZJ5j3HOSUlRRY7ffq05YkJiCNlIiKBsCkTEQmETZmISCBsykREAuFEXx3Qq1cvWUxpwuTYsWOyGB/YSbVh5MiRslhgYGCV95eWlladdITGkTIRkUDYlImIBMKmTEQkEDZlIiKBcKLPjtzptoZ+fn5mbf/999/LYoWFhdXKicgcY8aMser+UlNTrbo/kXCkTEQkEDZlIiKBsCkTEQmETZmISCCc6LMjMTExivH27dubtb3SRB+RNW3btk0x7u7uLouZe5vOJk2aVCclu8ORMhGRQNiUiYgEwqZMRCQQi5rygvnz0LljNHy9PREc6I8BTz2J06dOmaxTUlKChPgJCPDTwMdLjcEDB0Cn01k1aaI7YY2SvbNooi8z82s8P3482kW2w82bN/HGG6+hT59eOPrjT2jUqBEAIClxMr74Ig0fpmyBWq3GywnxGDTwaezbn1kjJ1CfPPHEE2avW1BQIIvt3bvXitmIiTVae3x8fGSxDh06KK6rNKln7kRffWNRU96ZanoP0/fWrENwoD+yso6gS5e/obCwEOvfX4cNGz9A166V3xR49721aNumNQ4ePIAOHaKtlzmRAtYo2btqXVPW/3HfBG8vbwBAVtYR3LhxAzGxscZ1Wt5/P4JDQnDwwAHFfZSWlkKv1xuXoiJ9dVIiMsEaJXtT5aZcUVGBxMRJ6NipEx784/HhulwdnJ2d4enpabKuVqOBLjdXcT8L5s+DxtfbuDQLC61qSkQmWKNkj6rclBPiJ+DEiRPY9MHmaiUwZeo05OVfNS5nc85Va39Ef2KNkj2q0i/6JibEIy1tF77am4GgoCBjXOunRVlZGQoKCkxGIrq8PGjvcHtJlUp1x1tSkqkuXbooxh0cHGSxrVu3ymLZ2dlWz0lUrNGaN27cOFmscePG1drnqlWrZLG8vLxq7dPeWDRSliQJExPisWP7NuzevQdhYWEm70dEROK+++5DRnq6MXb61Cn8dv48OkRzAoVqHmuU7J1FI+WE+An4aEsKPtn6Gdzc3ZH7xzU4tVoNV1dXqNVqjBg5ClOmJMLL2wseHh6YNDEB0dHRnNWmWsEaJXtnUVN+953VAIDu3WJN42vWYtiw4QCAhYsWw9HREf8YNBClpaWIi+uOpcnLrZQu0d2xRsneWdSUS8pu3nMdFxcXLF2WjKXLkqucFFFVsUbJ3vHWnYLq1KmTLBYaGqq4rtIvo86cOWPtlKgeCwgIkMVGjx5drX1mZWXJYpMnT5bFysrKqnUce8MbEhERCYRNmYhIIGzKREQCYVMmIhIIJ/oENWrUKFnM2dnZBpkQAa3/uHfI7Sx5dp5eL7+J05AhQ2Sx+japp4QjZSIigbApExEJhE2ZiEggbMpERALhRJ+g+vTpU63tDx8+bKVMqL5Ruk1pUlKSLGbJM/ZSUlJksdOnT1uWWD3BkTIRkUDYlImIBMKmTEQkEDZlIiKBcKJPUE8//bQstnPnTsV1582bJ4t9++23Vs+J6geNRiOLde3atVr7TEtLq9b29QlHykREAmFTJiISCJsyEZFA2JSJiATCiT5BKU3UeXt72yATqm+KiopkMaVnPjZv3lwWU3ruHgCkpqZWP7F6giNlIiKBsCkTEQmETZmISCBsykREAuFEHxGZKCgokMVatWpV+4nUUxwpExEJhE2ZiEggbMpERAIRrilb8ogZotvVVu2wRqmqzKkd4ZqywSD/NRGROWqrdlijVFXm1I5DSdlNof7Zr6iowKVLF+Hm5g6DoQjNwkJxNucc3N09bJ1alRUV6evEeQBinoskSTAYiuDvHwBHx5ofZ7BGxSbiuVhSo8J9Jc7R0RGBgUEAAAcHBwCAu7sHPDzE+HCro66cByDeuajV6lo7FmvUPoh2LubWqHCXL4iI6jM2ZSIigQjdlFUqFaa/9jpUKpWtU6mWunIeQN06F2uoK59HXTkPwP7PRbiJPiKi+kzokTIRUX3DpkxEJBA2ZSIigbApExEJhE2ZiEggwjbl1atWokXzZlC7N0KXzh1x+PAhW6d0T5mZX+PJvz+OsCbBcHFugB3bt5u8L0kSZs2cgdCQIHh6uKFXz+7IVnhKsK0tmD8PnTtGw9fbE8GB/hjw1JM4feqUyTolJSVIiJ+AAD8NfLzUGDxwAHQ6nY0ytg3WqO3U5RoVsil/8vHHmJKUiOmvvY4DBw/joTZt0a9Pb+Tl5dk6tbu6VlyMh9q0wdtLkxXfX7xoIVauWI7k5SuR+c13aNSwEfr27Y2SkpJazvTuMjO/xvPjx+PrzG+xK+3fuHHzBvr06YXi4mLjOkmJk7FrVyo+TNmCPXvTcenSRQwa+LQNs65drFHbqss1KuT3lLt07ojIdlF4e+kyAJU3gAlvGorxL7yIpClTbZydeVycG+DjT7ai/+OPA6gcgYQ1CUbCxJfx8qTJAIDCwkKEBAXgvTXrMHDQIFume1eXL19GcKA/9uxNR5cuf0NhYSGCAvywYeMHePKppwAAp37+GW3btMb+zG/QoUO0jTOueaxRsdSlGhVupFxWVoasrCzExMQaY46OjugaE4uDBw7YMLPqycnJQW5ursl5qdVqRLVvj4MHxT4vfWEhAMDbyxsAkJV1BDdu3EBM7K1zaXn//QgOCbHrPyNzsUbFU5dqVLimnJ+fj/Lycmi0GpO4VqOBTpdro6yq78/cNVqtSVyr0UKXK+55VVRUIDFxEjp26oQHW7cGAOhydXB2doanp6fJulqNRuhzsRbWqFjqWo0Kd+tOEktC/AScOHEC6Rn7bZ0KkaK6VqPCjZR9fX3h5OSEPJ3phIkuLw9arZ+Nsqq+P3PP+8vsry5PB62fmOc1MSEeaWm7sPvLrxAUFGSMa/20KCsrkz2KXpeXJ+y5WBNrVBx1sUaFa8rOzs6IiIhARka6MVZRUYF9GenoEC3uxfl7CQsLg5+fn8l56fV6HD50SLhJB0mSMDEhHju2b8Pu3XsQFhZm8n5ERCTuu+8+ZKTfOpfTp07ht/Pn7frPyFysUduryzUq5OWL+ISXMWb0SERERCIqKgrJyctQXFyMYcNH2Dq1uzIYDDibnW18fe5cDn48ehRe3t4ICQnBSxPiMW/uHISHN0doaChmzZwB/4AA4+y3KBLiJ+CjLSn4ZOtncHN3R+4f1+DUajVcXV2hVqsxYuQoTJmSCC9vL3h4eGDSxARER0cL95e3prBGbasu16iQX4kDgFUrV2DJksXQ5eaibdu2WPyvt9G+fQdbp3VX+/fvQ4+4brL4kKHDsGbtOkiShNmzZmLd2jUoKChAp86dsWzZcjRv0cIG2d6Zi7Pyv9XvrlmLYcOGA6j8Yv7UKUn4+KMtKC0tRVxcdyxNXg4/wf/X0JpYo7ZTl2tU2KZMRFQfCXdNmYioPmNTJiISCJsyEZFA2JSJiATCpkxEJBA2ZSIigbApExEJhE2ZiEggbMpERAJhUyYiEgibMhGRQP4fbbiHR3BFuwwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0399\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWUAAADOCAYAAADmKRMiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAa5ElEQVR4nO3deVhUZf8G8FtQBmQZNhkQRHF/raRAFDUrwX0pK5c2NZcWK8F9SSu3XEpNxSVNLZeULDV3zRIV316XVzTNNy39oZbKIBqrATqc3x/W2PQc9QAD88xwf65rrituzpnzHPj65TTPPGcq5RfeUkBERFJwsvUAiIjoDjZlIiKJsCkTEUmETZmISCJsykREEmFTJiKSCJsyEZFE2JSJiCTCpkxEJBE2ZSIiibApF8PKlSvg6lIZR4/+19ZDKVNLFn+MF57rhbp1wuDqUhkDB/S39ZBIo4pQo7/++iumTJ6ER1tEIzDAH8FBBrRtE4PvvvvW1kOzCjZlEsyc+SH27k1Co0aNULlyZVsPh8jC1i2bMWvmh6hdpy4mTJyEsW+PQ25OLjp37IAVKz6z9fBKjf/iSLD72z0IDQ1FpUqV4Oejt/VwiCw8/vgT+OVcKvz9/c3ZK6++hqZNIjFp4gT07fuy7QZnBbxSLqWBA/rDz0ePixcv4uluT8LPR4/atULx8aKFAIAfT55E+3Zt4OvthXp1ayNx7VqL/a9fv44xo0ci8pGH4eejRzU/HzzZtTNO/PCDcKwLFy7g2ae7wdfbCzWCgzByxDDs/mYXXF0qY9++vRbbHj58CF27dEKAvy989J5oE9sa33//b03nVLNmTVSqVKlkPxCSjqPVaKMHHrBoyACg0+nQoWMHXPrtN+Tk5BTzJyQXNmUrMJlMeKprF4SE1MDUadMRWrMmhsTHYeXKFXiya2dERkTi/anT4OnhiQH9X0Zqaqp539TU/8PmzZvRsVMnzPhwJoYOG45TP/6Itm1icPnyZfN2eXl56NC+Lfbs+Q5vvPkWRo8Zi4P/+Q/Gvf22MJ6kpD1oE9MaOdnZGDf+HUyaPAVZmVno0K4tjhw5XC4/E5JLRajRtDQjqlatiqpVq5Zof2nkF95S+ND2WLJ0mQJA+fd/Dpqzl3r3UQAokyZPMWdp6RmKm5ubUqlSJWXV6jXm/MTJUwoAZdz4d8xZVk6eciO/0OI4p38+q+h0OuXd9yaYsxkffKgAUL78aoM5y8zOVRo0aKgAUHbt/lbJL7yl/FFwU6lbt57Stl075Y+Cm+Ztf8/KUWqFhSmxbdoU65zd3d2Vl3r3sfnPng/W6L0ep/53WnF1dVVeePElm/8OSvvglbKV9Os/wPzf3t7eqF+/Adzd3dG9Rw9zXr9BA3h7e1tcheh0Ojg53f41mEwmXLt2DR4eHqhfvwGOHztm3u6bb3YhODgYXbp2NWeurq7oP+DOcQHgh+PHcfbsL+jV63lcu3YNGRkZyMjIQF5eHlq3jsGB5GQUFRVZ/fxJfo5aozdu3MALzz8HNzc3THl/qvYfiKQ40WcFrq6uqFatmkXmpfdCcHCI8Nqsl16PzN9/N39dVFSE+QnzsHjxxzifmgqTyWT+nq+fr/m/L164iLDatYXnq1OnrsXXZ8+eBQAMHNDvruPNysqCj4+PxrMjR+CoNWoymdD7xRfw00//w6YtW1G9evX77iM7NmUrcHZ2Llau4M4ncM2YPg0TJ7yHvi/3w3vvTYSPrw+cnJwwcvjwEl3R/rXPtOkz0Dg8XHUbDw+PYj8v2TdHrdFBr7+G7du34bMVq9C6dUyxxyIjNmUb27hhAx5/4gksXvKJRZ6VlQk/fz/z16E1Q3H6p5+gKIrFlci5c2ct9qtdpzYAwNPLC7Gxbcpw5FRRyFqjY8eMwsoVn2HmrNno9dxzJX4e2fA1ZRtzdnaGolh+du36r77CpUuXLLK2bdvh0qVL2LpliznLz8/H8mXLLLaLiIhE7Tp1MGf2bOTm5grHu3r1qhVHTxWBjDU6e9ZMfDR7NkaNHoO3BscV53SkxytlG+vYqROmvj8FrwwcgOjmzXHqxx+RuHYNwmrXtthu4CuvYtGihejT+0W8+dZgBAUFIXHtGri6ugKA+crEyckJiz5ejKe6dsEjDzdGnz59UT04GJcvXcK+fXvh5emFDV9vuueYtm3dghMnTgAAbt68iR9PnsS0qe8DALp06YqHGje29o+BJCZbjW76+mu8PXYM6tath4YN/4U1n39u8f3YNm1gMBis/FMoP2zKNjZ6zFjcuJGHLxIT8dWX6/DwI49g46bNGD/O8r2dHh4e2LlrN4YNjceC+Qnw8PDAiy+9hOjoFniuVw9z4QO3Vzzt238A06a+j48XLURubi4MgYFoGtUUA1555b5j2rhxI1avWmn++vjxYzh+/PYse3BICJtyBSNbjZ44cXvRytmzv6B/v77C93ft/taum3Kl/MJbyv03I1klzJuLkSOG41zqBQQHB9t6OEQC1mjxsCnbkT/++ANubm7mr/Pz89EsqglMJhN+/N9PNhwZ0W2s0dLjyxd2pFfP7qhRIxTh4eHIysrC2jVrcObMaXy2YpWth0YEgDVqDWzKdqRt23b4dPlyJK5dA5PJhH/9qxFWrV6DHj172npoRABYo9bAly+IiCTC9ykTEUmkzF6++HjRQsyePQvGtDQ0btwYs+fMRVRU0/vuV1RUhCtXLsPDw5P39CVNFEVBbm4OgoKqm2+cowVrlMpLsWq0LG49t2r1GsXFxUVZ/MlS5djxE0r/AQMVb29v5eJvl++777nU8woAPvgo9uNc6nnWKB9SP7TUaJm8ptyqZXNENonCnLnzANy+sqhbuxYGvfEmRo4afc99s7KyYKjmd89tiNQYr16DXq/XtC1rlGxBS41a/TXlwsJCpKSkICYm9s5BnJzQOiYWhw4eFLYvKChAdna2+ZGba98f5UK2o/WlBNYo2YqWGrV6U87IyIDJZEKAIcAiNwQEwGhME7b/YMZ0BPj7mh91wmpZe0hEFlijJDObv/ti1OgxSM+4bn6cSz1v6yERWWCNUnmy+rsv/P394ezsjHRjukVuTE+HwRAobK/T6aDT6aw9DKK7Yo2SzKx+pezi4oKIiAgkJe0xZ0VFRdibtAfNoqOtfTiiYmONkszK5H3KcfFDMXBAP0RERCIqKgoJCfOQl5eHPn1fLovDERUba5RkVSZNuUfPnsjIuIpJkybAmJaG8PBwbN66za7vcUqOhTVKspLu3hfZ2dkI8Pe9/4ZE/5CecR1eXl5lfhzWKJWUlhq1+bsviIjoDjZlIiKJsCkTEUmETZmISCJsykREEmFTJiKSCJsyEZFE2JSJiCTCpkxEJBE2ZSIiibApExFJhE2ZiEgibMpERBJhUyYikgibMhGRRMrkJvdUPE5O4t/GKlWqaN7/hRdeELKJEycKWUhIiJA1adJEyE6dOiVkb7zxhpD5+or3FB4yZIiQLV68WMj++OMPIZswYYKQmUwmISNyZLxSJiKSCJsyEZFE2JSJiCTCpkxEJBF+cGoZcnV1FTK1T0t+9913haxfv36aj6M2GVZYWChkV65cEbKTJ08K2VNPPaX52NYUGxsrZElJSZr35wenkuz4walERHaGTZmISCJsykREEmFTJiKSCFf0WUmPHj2ETG0C74EHHrD6sdWO8/nnnwvZyJEjhez111/XdIzMzEwhmzNnjpDVq1dPyLp16yZk7u7uQubv769pLGTpscceEzK1VZ5qXn31VSHbsGGDkGVkZAjZ6dOnhWznzp2qx1HbltTxSpmISCJsykREEmFTJiKSCJsyEZFEuKKvBAIDA4Vs9+7dQqZ1Uu/WrVtCdvXqVSFLTExU3f+jjz4Sst9++03Tsd955x0hU5s46t69u5BlZWVpOsY333wjZG3atNG0r9ptTe+moq7o27Fjh5C1a9dOyBRF/KdeqVIlq26ndktWANi4caOQ9enTR3VbR8YVfUREdoZNmYhIImzKREQSYVMmIpIIV/Tdh06nEzK1FU9aJ/WWLVsmZNu3bxcytYmRsjB58uQS7+vp6anp+dQmDtVoXV1YkdWsWVPIIiIihExtYk6NtbdTW6kJAC+++KKQdejQQcg6duwoZEePHtV0bEfBK2UiIomwKRMRSYRNmYhIImzKREQSKfZEX3Lyfnw0axaOHUvBlStXsO7L9Xjyb5/ppigKJk2cgE+XL0NmZiaat2iBhIQFqKtyS0d7UFBQIGSHDx8WsujoaE3Pt3//fiErr0m90lBbhbRu3TohU1tJpmbVqlVCtnTp0uIPTIUj1+jAgQOFzM/PT8jUVoT27t1byNRuyammWrVqQjZ27Fgha9Wqler+aqv/1Ma9bds2IXviiSeEzJFvBVrsK+UbeXl4qHFjzJmboPr9WTM/xMIF85EwfyGSD3wP96ru6NKlE/Lz80s9WCItWKNkz4p9pdy+Q0e07yC+bQW4/ddwfsI8jBn7Nro++SQAYNmnnyE0pDo2b9qEnr16CfsUFBRYXI3m5GQXd0hEFlijZM+s+ppyamoq0tLSEBNz56Pi9Xo9opo2xaFDB1X3+WDGdAT4+5ofdcJqWXNIRBZYoyQ7qzZlozENABBgMFjkhgADjGlpqvuMGj0G6RnXzY9zqeetOSQiC6xRkp3NV/TpdDrVVXOycHZ2FjK1CQp7oLYCT22FmNrqxKFDhwpZ06ZNhezMmTNCNn36dCFT+wzBoqIiIZOBTDU6btw4IVObRFO7XapaVhq7du0Ssn379qlu++ijjwqZ2ipBtQlFtRW0aiv/Lly4oHpse2PVK2WD4fZ9htONRovcmG6EQeUexETljTVKsrNqUw4LC0NgYCCSkvaYs+zsbBw5fBjNmml7yxhRWWKNkuyK/fJFbm4uzp09a/76/PlU/HD8OHx8fREaGoq3Bsdh+rSpqFu3HmrVqoWJE95DUPXqFu8TJSpLrFGyZ8VuykeP/hft2975KJ9RI0cAAF7q3QdLly3H8BEjkZeXhzffeB2ZmZlo0bIltmzZBldXV+uNmugeWKNkz/gZffdRubL4d+v5558XshUrVmh6vq1btwrZk3++X9aa6tevL2R16tQRMrVbiap9BuH169eFTG1ipX379kKmddVYaVWEz+gzmUxCpjbRd+zYMSGLiooqkzH9ndrvH1BfqVeazwdcsmSJkA0aNEjLEG2Kn9FHRGRn2JSJiCTCpkxEJBE2ZSIiiXCirwRcXFyEbP369ULWuXNnIcvLy9O0ndotPu+mUaNGQqa2Yi48PFzzc/6T2qq8t99+u8TPVxYqwkRfenq6kKmtML127ZqQTZ06Vch27typad+qVasKmb+/v5Cp3ZIVABo0aCBkpZnoU9uudevWQlacf0flgRN9RER2hk2ZiEgibMpERBJhUyYikojNb91pjwoLC4XswIEDQqY2gefu7i5kMTExQhYUFCRkap/PBgCPPPKIkPn6apuISklJEbJDhw4J2ZQpUzQ9H5WtYcOGCdnMmTOFTO0WmGrbTZ48WcjUVmCqTfSpTTCqTcoB6hNzaqvynn76aSFTOxe15+vWrZuQyTbRpwWvlImIJMKmTEQkETZlIiKJsCkTEUmETZmISCJcZm0lakuv1Wa733rrrfIYjiq1pddq96DNzc0tj+FYXUVYZq2mYcOGQqb2YaPlsdT5bu++UBtP9+7dhezIkSNCFhkZqenYau8kKo97SBcHl1kTEdkZNmUiIomwKRMRSYRNmYhIIlxmbSVOTuLfNw8Pj3I5ttqyaLUJE7X7H9vrpB7dcfr0aSFr0qSJkKlNCD722GNCpjYhqNXdPiR32rRpJX5OtUk9tcxR8EqZiEgibMpERBJhUyYikgibMhGRRDjRZyUPPvigkD388MNWPcb169dV89dee03ITpw4YdVjk325ceOGkKmteFPLbCk5OVnI1Fb0OTJeKRMRSYRNmYhIImzKREQSYVMmIpIIJ/qspF69ekKmNtGntvru999/F7IOHToImZubm+qx27dvL2Sc6CN7pLY6kSv6iIjIZtiUiYgkwqZMRCQRNmUiIolwou8+qlSpImTLly8Xsi5dugjZpUuXhExtUm7GjBmaxnLmzBnVfMGCBZr2J7JHd/vcv39SW/kXEREhZLKtYvwnXikTEUmETZmISCJsykREEilWU/5gxnS0bB4Nf19v1AgOQo9nn8HP/3idMz8/H/Fxg1E9MAB+Pno817MHjEajVQdNdDesUbJ3xZroS07ej9cGDUKTyCa4desW3n13PDp37ojjP5yEu7s7AGDkiOHYsWM7Pl+bCL1ej6HxcejVszv27hNvyScbZ2dnIVu9erWQ9ejRQ9PzTZo0Sciys7OF7MKFC5qer1GjRqq52mTGgQMHND2no3H0Gq2IKtqKvmI15S1bt1t8/cnS5agRHISUlKNo1eoxZGVl4bNPl2PFytVo3ToGALDkk2UIb/wgDh06iGbNoq03ciIVrFGyd6V6TTk7KwsA4OvjCwBISTmKmzdvIiY21rxNg4YNUSM0FIcOHlR9joKCAmRnZ5sfOTnilSRRSbFGyd6UuCkXFRVhxIhhaN6iBR7481M3jGlGuLi4wNvb22JbQ0AAjGlpqs/zwYzpCPD3NT/qhNUq6ZCILLBGyR6VuCnHxw3GqVOnsGr1mlINYNToMUjPuG5+nEs9X6rnI/oLa5TsUYlW9A2Jj8P27dvw7XdJCAkJMeeGQAMKCwuRmZlpcSViTE+HITBQ9bl0Oh10Ol1JhmF1iYmJQvbss89q2jc9PV3IFi9erGnftLtcof2T2i0+ASAjI0PT/hWJo9ZoRaR1RZ/adtWqVbP2cMpcsa6UFUXBkPg4bN70NXbt2o2wsDCL70dERKJKlSpI2rPHnP185gx+vXgRzaI5gUJljzVK9q5YV8rxcYPxReJafLl+Azw8Pc1XeHq9Hm5ubtDr9Xi5X3+MGjUCPr4+8PLywrAh8YiOjuasNpUL1ijZu2I15SWLPwYAtGsTa5kvXYY+ffoCAD6cOQtOTk54vldPFBQUoG3bdpibMN9KwyW6N9Yo2btiNeX8wlv33cbV1RVz5yVg7ryEEg+KqKRYo2TveOvOv+nevbuQaV05pDah8Ouvv2ra96+VZvfj5+enmquNe8qUKZqek0h2pVnR17BhQyHbtWtXqcdUlnhDIiIiibApExFJhE2ZiEgibMpERBLhRN/fdO3aVcjGjx8vZE2bNhUytdVEPj4+1hnYnypXVv91Pf7440LGiT5yFFpX9Dk5ideYycn2dztWXikTEUmETZmISCJsykREEmFTJiKSCCf6/mbr1q1CtmPHDiFTm1Do1KmTkEVGRlpnYPexcuXKcjkOkS1oXdFXVFRUHsMpc7xSJiKSCJsyEZFE2JSJiCTCpkxEJBFO9N2HyWTSlG3atElTRkTFo3VFX0pKiqZMdrxSJiKSCJsyEZFE2JSJiCTCpkxEJBFO9BGRNDZs2CBkNWrUELJu3bpp2tce8UqZiEgibMpERBJhUyYikgibMhGRRCrlF94S74FnQ9nZ2Qjw97X1MMgOpWdch5eXV5kfhzVKJaWlRnmlTEQkETZlIiKJsCkTEUlEuqas9jEvRFqUV+2wRqmktNSOdE05NzfH1kMgO1VetcMapZLSUjvSvfuiqKgIV65choeHJ3Jzc1AnrBbOpZ6Hp2fZz6qXlZycbIc4D0DOc1EUBbm5OQgKqq76obbWxhqVm4znUpwale7eF05OTggODgFw5+bWnp5e5fJWp7LmKOcByHcuer2+3I7FGrUPsp2L1hqV7uULIqKKjE2ZiEgiUjdlnU6HcePfgU6ns/VQSsVRzgNwrHOxBkf5eTjKeQD2fy7STfQREVVkUl8pExFVNGzKREQSYVMmIpIImzIRkUTYlImIJCJtU/540ULUr1cHek93tGrZHEeOHLb1kO4rOXk/nun2FMJq1oCrS2Vs3rTJ4vuKomDihPdQKzQE3l4e6NihHc7+8ouNRnt3H8yYjpbNo+Hv640awUHo8ewz+PnMGYtt8vPzER83GNUDA+Dno8dzPXvAaDTaaMS2wRq1HUeuUSmb8pfr1mHUyBEYN/4dHDx0BA81DkfXzp2Qnp5u66Hd0428PDzUuDHmzE1Q/f6smR9i4YL5SJi/EMkHvod7VXd06dIJ+fn55TzSe0tO3o/XBg3C/uR/Y9v2nbh56yY6d+6IvLw88zYjRwzHtm1b8fnaROz+bg+uXLmMXj2723DU5Ys1aluOXKNSvk+5VcvmiGwShTlz5wG4fQOYurVrYdAbb2LkqNE2Hp02ri6Vse7L9XjyqacA3L4CCatZA/FDhmLosOEAgKysLISGVMcnS5ejZ69ethzuPV29ehU1goOw+7s9aNXqMWRlZSGkeiBWrFyNZ559FgBw5vRphDd+EPuSD6BZs2gbj7jssUbl4kg1Kt2VcmFhIVJSUhATE2vOnJyc0DomFocOHrThyEonNTUVaWlpFuel1+sR1bQpDh2S+7yys7IAAL4+tz+XLiXlKG7evImY2Dvn0qBhQ9QIDbXr35FWrFH5OFKNSteUMzIyYDKZEGAIsMgNAQEwGtNsNKrS+2vsAQaDRW4IMMCYJu95FRUVYcSIYWjeogUeePBBAIAxzQgXFxd4e3tbbGsICJD6XKyFNSoXR6tR6W7dSXKJjxuMU6dOYU/SPlsPhUiVo9WodFfK/v7+cHZ2RrrRcsLEmJ4OgyHQRqMqvb/Gnv6P2V9juhGGQDnPa0h8HLZv34Zd33yLkJAQc24INKCwsBCZmZkW2xvT06U9F2tijcrDEWtUuqbs4uKCiIgIJCXtMWdFRUXYm7QHzaLlfXH+fsLCwhAYGGhxXtnZ2Thy+LB0kw6KomBIfBw2b/oau3btRlhYmMX3IyIiUaVKFSTtuXMuP585g18vXrTr35FWrFHbc+QalfLli7j4oRg4oB8iIiIRFRWFhIR5yMvLQ5++L9t6aPeUm5uLc2fPmr8+fz4VPxw/Dh9fX4SGhuKtwXGYPm0q6tath1q1amHihPcQVL26efZbFvFxg/FF4lp8uX4DPDw9kfbna3B6vR5ubm7Q6/V4uV9/jBo1Aj6+PvDy8sKwIfGIjo6W7h9vWWGN2pYj16iUb4kDgEULF2D27FkwpqUhPDwcsz6ag6ZNm9l6WPe0b99etG/bRshf6t0HS5cth6IomDRxApYvW4rMzEy0aNkS8+bNR7369W0w2rtzdVH/W71k6TL06dMXwO035o8eNRLrvkhEQUEB2rZth7kJ8xEo+f8aWhNr1HYcuUalbcpERBWRdK8pExFVZGzKREQSYVMmIpIImzIRkUTYlImIJMKmTEQkETZlIiKJsCkTEUmETZmISCJsykREEmFTJiKSyP8DV9ljWv4piLgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "test_dataloader = DataLoader(siamese_test, shuffle=True, num_workers=8, batch_size=1)\n",
    "\n",
    "def show_image_pair(imgA, imgB, label, similarity_score, i):\n",
    "    \n",
    "    fig, ax = plt.subplots(1, 2, figsize=(4, 4))\n",
    "    ax[0].imshow(imgA.squeeze(), cmap='gray')\n",
    "    ax[0].set_title('Image 1')\n",
    "    ax[1].imshow(imgB.squeeze(), cmap='gray')\n",
    "    ax[1].set_title('Image 2')\n",
    "    plt.savefig(f\"image_{i}.jpeg\", bbox_inches=\"tight\", dpi = 300)\n",
    "    print(similarity_score)\n",
    "    plt.show()\n",
    "\n",
    "def visualize_siamese_pairs(data_loader, total_images=4):\n",
    "    for idx, batch in enumerate(data_loader):\n",
    "\n",
    "        if idx == total_images:\n",
    "            return\n",
    "            \n",
    "        imgA, imgB, label = batch\n",
    "        outputA, outputB = net(imgA.cuda(), imgB.cuda())\n",
    "        euclidean_distance = F.pairwise_distance(outputA, outputB)\n",
    "        similarity_score = torch.exp(-euclidean_distance)\n",
    "        \n",
    "        imgA = imgA[0].numpy()\n",
    "        imgB = imgB[0].numpy()\n",
    "        label = label[0].item()\n",
    "        \n",
    "        show_image_pair(imgA, imgB, label, round(similarity_score.item(), 4), idx)\n",
    "        \n",
    "visualize_siamese_pairs(test_dataloader, total_images=4)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
